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ABSTRACT 

We report the detection of far-IR CO rotational emission from the prototypical Seyfert 2 galaxy NGC 
1068. Using Herschel-PACS, we have detected 10 transitions in the Juppor = 14 — 24 (i?uppcr/fcs — 
580 — 1656 K) range, all of which are consistent with arising from within the central 10" (700 pc). 
The detected transitions are modeled as arising from 2 different components: a moderate excitation 
(ME) component close to the galaxy systemic velocity, and a high excitation (HE) component that is 
blueshifted by ^70 kms~*^. We employ a large velocity gradient (LVG) model and derive nH2 ^ 10^ '' 
cm-^ Tkin - 150 K, and Mh2 ^ lO^'^ Mq for the ME component, and nH2 ^ lO*^'^ cm^^ Tan ^ 440 
K, and AIh2 ~ 10^'^ Mq for the HE component, although for both components the uncertainties 
in the density and mass are ~±1 dex. Both components arise from denser and possibly warmer gas 
than traced by low- J CO transitions, and the ME component likely makes a significant contribution 
to the mass budget in the nuclear region. We compare the CO line profiles with those of other 
molecular tracers observed at higher spatial and spectral resolution, and find that the ME transitions 
are consistent with these lines arising in the ~200 pc diameter ring of material traced by H2 1-0 S(l) 
observations. The blueshift of the HE lines may also be consistent with the bluest regions of this 
H2 ring, but a better kinematic match is found with a clump of infalling gas ~40 pc north of the 
AGN. We consider potential heating mechanisms, and conclude that X-ray or shock heating of both 
components is viable, while far-UV heating is unlikely. We also report sensitive upper limits extending 
up to Juppor — 50, which place constraints on the emission from the X-ray obscuring medium. 
Subject headings: galaxies: active — galaxies: individual(NGC 1068) — galaxies: ISM — galaxies: 
nuclei — galaxies: Seyfert — infrared: galaxies 



1. INTRODUCTION 

The excited molecular gas in the centers of Seyfert 
galaxies offers a sensitive probe of the nature of ac- 
tive galactic nucleus (AGN) feedback on the surround- 
ing ISM. Observational studies of this materi al in Seyfert 
nucle i have typically used the H2 rotational ( Lutz et al.l 
[2OOOI: iRigopoulou etaLl I200I IRoussel et al.l I2007D and 
the well-studied ro- vibrational (e. g., [T hompson et al.l 

[ lDavi( 
transitions. The pure H2 
> 500 K) are easily ther- 
densities, while 



the ro-vibrational lines (-Euppor/^B > 7000 K) may be 
excited through collisions in dense (nH2 > 10^ cm" 



19781: .Mouri 1991 iMaloneyl 119971: iDavic s ct al. 2005,; 
Rodriguez- Ardila et all 120051 ) 



rotational lines (-Euppor/^B 
malized at moderate (fiH2 
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gas o r through UV fiuorescence ([Sternberg fc Dalgarnol 
119891) . Observations of these tracers in Seyferts have 
identified a number of potentially important e xcitation 
mechanisms, including X-rays from the AGN (jMalonevI 
[19971 : IRodn'g uez- Ardila et al.l 120051 ). shocks associated 
with superno va remnants, r adio j ets, and gravitational 
instabilities (IRoussel et al.l l2007f). and stellar far-UV 
(FUV) radiation pavies et al.ll2'005[ ). with no clear con- 
sensus on a single dominant excitation source. The 
far-IR (FIR) CO rotational transitions (CO[Juppcr — > 
Juppor — 1], with Juppor ~ 13 — 50) arisc from states 
500 — 7, 000 K above ground and have critical densities 
of ^ 10^ — 10^ cm~^, and complement the H2 transi- 
tions for studies of warm and dense material. Com- 
pared with H2, the FIR CO lines trace similar energy 
levels, but have higher critical densities, and are less 
sensitive to extinction. Additionally, the smaller en- 
ergy gaps between levels leads to a finer sampling of 
density-temperature phase space. These lines have been 
proposed as potential too ls for studying the o bscuring 
medium of type 2 AGN ([Krolik fc LeppI 119890 . deter- 
mining t he energy budgets of c omposite starburst/AGN 
systems (iMeijerink et al.l 120071 ). and identifying accret- 
ing b l ack holes in the early u niverse (^paans & Mciicrin3 
[20081 ISchleicher et al.ll2010l) . but previous facilities were 
unable to detect this line emission from extragalactic 
sources. Here we take advantage of the superb sensitiv- 



2 



Hailey-Dunsheath et al. 



ity of Herschel-PACS to conduct the first extragalactic 
study of FIR CO emission, from the prototypical Seyfert 
2 galaxy NGC 1068. 

NGC 1068 is one of the brightest and best stud- 
ied Seyfert 2 galaxies. The paradigm of an opti- 
cally and geometrically thick molecular torus account- 
ing for the Seyfert type 1 and 2 dichotomy followed 
the det ection of scattered broad l ine emission from this 
source (IMiller fc Antonuccil [l98l . and NGC 1068 has 
been at the center of subsequent studies of the ISM in 
Seyfert nuclei. The molecular gas in the central 1' of 
NGC 1068 has been well studied, and here we review 
some of the key results. Interferometric observations 
of CO (1-0) have id entified a pair of » 1 5" (f^l.l kpc) 
radius spi ral arms (iPlanesas "etaII[T99l IHelfer fc Bhtl 
[1995; Sch innerer et al J12OOO0, which m ay be modeled as 
forming in response to a sal 7 kpc bar (" Schinnerer et aF 



200C). These arms are bright in 61-7 ([Davies et al 



T998I) . PAH emission (iLe Floc'h eFall I2001D. an d sub- 



millimeter continuum ()Papadopoulos fc Seaauistlll999aD . 
and contain most of the star formation in the cen- 
tral region. Centered on the AGN is the ~5" 
(^350 pc) circumnuclear disk (CND), which is visi- 
ble in CO and H2 1-0 S(l), bu t becomes particu- 
larly prominent in images of HCN (iTacconi et al.l[T99l 
and other high density tracers ()Garcfa-Burillo et al.l 
IMp). Strong emission in CO(4-3) and HCN(1- 
0) ind icate the gas in the CND is both warm and 
dense (jTacconi et al.1 11994': 'Sternberg et al.l 119941 : [Israeli 
120091 ) (although see iKri ps et al.. (20 111 ) for a lower 
density model). The high abundances of HCN, 
CN, HsO"*", and other molecules in the CND sug- 
gest an X-ray— driven chemistry (lUsero et al.l 120041 : 
iGarcfa-Burillo et al.llMol : lAalto et al.ll2011[ ). and X-ray 
heating has also been invo ked to explain the strong H2 1- 
S(l ) and [Fe II] emissio n (iRotaciuc et al.lll991l : lMalonevl 
[19971 IGalliano fc Alloinll2002D . At ^0.3" resolution the 
H2 1-0 S(l) observ ations resolve the CND into a ~1" 
ring-like structure ([Galliano fc AlloinI l2002f) , while the 
line spectral profiles show evidence fo r rotation, expan- 
sion, and more complex k inematics (Galliano fc Alloin 
[200a IGalliano et al.l l2003t iDavies et al. .2008, ). Shocks 
following this non-circular motion, and possibly associ- 
ated with jet-ISM interactio ns, may also be im portant in 
heating the molecular CND (iKrips et al.ll201lD . At -0.1" 
resolution the H2 1-0 S(l) images reveal two clumps of in- 
falling molecular material at ~0. 1 — 0.4" scales that likely 
play an important role in both fueling and obscuring the 
AGN, with an estimated infall ra te of —15 Mf^iyr"^ to 
withi n a few parsecs of the nucleus ([Miiller Sanchez et al.l 
[2OO9I) . Finally, milli -arcsec resolution radio observations 
identify a series of H2O maser spots that trace out the 
inner surface of a —0.5 pc radius molec ular disk, cen- 
tered on the AGN ([Gallimore et al.ll2004 and references 
therein) . 

Hard X-ray observations of NGC 1068 indicate 
large obscuration to the nucleus (llwasawa et al.l 119971 : 
IMatt et all [19971 : IColbert et all I2002t ) bv an interven- 
ing medium with a column density p ossibly exceeding 
A^H > 10^^ cm-2 ([Matt et al.lll997| ). Interferometric 
mid-IR observations of NGC 1068 identify a parsec scale 
structure of hot dust that, along with the H2O maser 
disk, may represent the dusty molecular torus resp onsible 
for the X-ray obscuration (iJaffc et al..,2004; .Raban et al.l 



|2009[) . However, other investigators have found evidence 
that at least some of the nuclear obscur ation occurs on 
few to ten parse c scales fro m the AGN dC ameron et al.l 
119931 iHonig et a l. 2006; Mii ller Sanchez e t al. 200^. 

The organization of this paper is as follows. In sec- 
tion [2| we describe the Herschel-PACS observations of 
the FIR CO lines in NGC 1068. In section O we an- 
alyze the gas excitation and estimate physical parame- 
ters, and in section |4| we compare the physical parame- 
ters and line profiles with those of other molecular gas 
tracers. In sections [5| and [6| we discuss potential heat- 
ing mechanisms, in section [7| we discuss non-detections, 
and in section [5| we summarize our findings. Through- 
out t his paper we adopt a distance to NGC 1068 of 14.4 
Mpc ([Bland- Hawthorn et al.l [19971) . and a systemic ve- 
locity t^LSR = 1125 kms~^. 

2. OBSERVATIONS 

The observations were made with the Photode- 
tector Array Cam era and Spectrometer (PACS; 
iPoglitsch et al.ll2010D on board the Herschel Space Ob- 
servatory (IPilbratt et an[201fl[ ). as part of the SHINING 
guaranteed time key program. Ten high resolution range 
scans were concatenated to cover the 52 — 98 /im and 
104 — 196 /im ranges, and we additionally obtained 
deeper integrations of CO(17-16), 00(24-23), and 
CO (40-39). These observations amounted to a total of 
11.5 hours of integration time. The data reduction was 
done using the standard PACS reduction and calibration 
pipeline (ipipe) included in HIPE 5.0 975. However, 
for the final calibration we normalized the spectra to 
the telescope fiux and recalibrated it with a reference 
telescope spectrum obtained from dedicated Neptune 
observations. 

In Figure [1] we show the spectra from the central spa- 
tial pixel (spaxel) centered on 11 of the 12 CO transitions 
falling in the 104-196 /im range. The 00(25-24) line at 
Aiest — 104.44 /im lies in a noisy region at the edge of this 
range, and is not included. The spatial distribution over 
the PACS array of each detected CO line is consistent 
with that of an unresolved source. All fluxes and upper 
limits presented here were therefore extracted from the 
central spaxel, and referenced to a point source by divid- 
ing by the recommended beam size correction factors. 
The CO line fluxes were measured by fitting a Gaussian 
profile to the baseline-subtracted spectra. In most cases 
the three free parameters of the Gaussian were allowed 
to vary, while baseline uncertainties necessitated fixing 
the widths of the 00(15-14) and 00(22-21) lines to typ- 
ical values derived from other line fits (corresponding to 
a deconvolved FWHM of 250 kms"^ see Figure[2|). The 
00(16-15) line is blended with the 163 /im OH doublet, 
and 00(17-16) with a pair of flanking 0H+ lines. In 
both cases we estimate the CO flux by simultaneously 
fltting all features. 00(23-22) is blended with a strong 
H2O line with a rest wavelength 209 kms~^ to the red. 
If the combined feature were attributed solely to H2 O it 
would be both broader and more blueshifted than any 
of the other 6 H2O lines detected in the PACS scans, 
and we interpret this as evidence for signiflcant contam- 
ination by 00(23-22). A comparison with the average 
of the CO(22-21) and CO(24-23) profiles suggests the 
00(23-22) line is weaker and/or more redshifted than 
these lines (Figure [2]). However, due to the uncertain- 
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TABLE 1 
PACS CO Line Observations 



Line 




Flux^ 


Vo 


1^ WHM 






[10 ' Wm ^] 


n — 1 1 
[kms 


n — 1 1 
[kms -^J 


CO(14-13) 


186.00 


7.2 ± 2.3 


19 ± 23 


305 ± 32 


CO{15-14) 


173.63 


6.8 ± 2.1 


2 ± 25 


329'' 


CO(16-15) 


162.81 


8.1 ± 2.5 


49 ± 25 


369 ± 28 


CO{17-16) 


153.27 


5.8 ± 1.8 


± 26 


359 ± 27 


CO(18-17) 


144.78 


5.1 ± 1.6 


— 15 ± 28 


379 ± 31 


CO(19-18) 


137.20 


2.6 ± 0.9 


-31 ± 35 


324 ± 55 


CO{20-19) 


130.37 


2.5 ± 0.9 


-62 ± 36 


297 ± 55 


CO{21-20) 


124.19 


2.4 ± 0.9 


-17 ± 46 


407 ± 87 


CO{22-21) 


118.58 


4.0 ± 1.4 


-44 ± 43 


387"= 


CO(23-22) 


113.46 


blended 






CO{24-23) 


108.76 


2.6 ± 1.0 


-94 ± 49 


385 ± 95 


CO{25-24) 


104.44 


< 11.2 






CO(26-25) 


100.46 


n/a 






CO{27-26) 


96.77 


< 5.8 






CO{28-27) 


93.35 


< 6.0 






CO(29-28) 


90.16 


< 9.3 






CO{30-29) 


87.19 


< 6.2 






CO{31-30) 


84.41 


blended 






CO{32-31) 


81.81 


< 7.4 






CO{33-32) 


79.36 


< 9.5 






CO{34-33) 


77.06 


< 5.8 






CO{35-34) 


74.89 


< 6.2 






CO{36-35) 


72.84 


< 8.3 






CO{37-36) 


70.91 


< 7.9 






CO{38-37) 


69.07 


< 10.1 






CO{39-38) 


67.34 


< 19.2 






r^r\( /in '3^^ 
CO(4U-o9j 


65.69 


< 22.3 






CO(41-40) 


64.12 


< 13.3 






CO{42-41) 


62.62 


< 21.7 






CO{43-42) 


61.20 


< 16.9 






CO{44-43) 


59.84 


< 10.8 






CO{45-44) 


58.55 


< 14.9 






CO{46-45) 


57.31 


blended 






CO{47-46) 


56.12 


< 14.8 






CO{48-47) 


54.99 


< 28.4 






CO(49-48) 


53.90 


< 31.9 






CO{50-49) 


52.85 


< 45.4 







"Total uncertainties combine a 30% calibration error with statis- 
tical errors in line fits. Upper limits are 3(T, and refer to the flux 
density integrated over a 600 kms~^ bin. Some lines are blended 
with a strong feature, and it is not possible to obtain a flux mea- 
surement or useful upper limit. CO(26-25) was not covered in the 
PACS scans. 

''Relative to Vlsr = 1125 kms~^. Total uncertainties combine 
a spectral calibration accuracy of 10% of the spectral resolution 
(20 — 30 kms~^) with statistical errors in line fits. 

"^Fixed to an intrinsic FWHM of 250 kms~^ (see text). 

ties involved in deconvolving the CO and H2O lines, we 
simply exclude 00(23-22) from our analysis. All other 
transitions from 00(14-13) through 00(24-23) are well 
detected (Table [1|). 

The 52 - 98 range includes the 00(27-26) through 
00(50-49) transitions, and none of these lines are de- 
tected. We estimate upper limits by first binning the 
data to 600 kms~^ bins, and then estimating the 3cr 
noise levels (Table [IJ . We also detect no emission from 
1300. The i3C0(14-13) transition at Arcst = 194.55 jim. 
lies at the noisy edge of a scan, while the Jupper > 21 
transitions at A > 104 /xm are blended with ^^00 lines. 
For 1300(15-14) through 1^00(20-19) we estimate Sc- 
upper limits of (2 — 4) X lO^^'' Wm~^, also for a 600 
kms^i bin size. 

3. EXCITATION ANALYSIS 
3.1. Evidence for Two Components 
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Fig. 1. — Continuum-subtracted spectra of the Jupper = 14 — 24 
CO lines. The blue curves show the line fits, and features other 
than CO are labeled in red. For CO(16-15) and CO(17-16) the 
green and red curves show the decomposition of the fit into CO and 
other lines, respectively. All lines are detected with the exception 
of CO(23-22), which is blended with a strong H2O line 209 kms"! 
to the red. Here the overplotted green curve is an average of the 
CO(22-21) and CO(24-23) profiles as a reference. 
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Fig. 2. — Top: FIR CO line fluxes and upper limits measured 
here, along with lower-J lines from the literature. The bright- 
est set of Juppcr = 1 — 4 line fluxes are measured in 11" — 21" 
beams that contain a mixture of the CND and the more extended 



emission (Israel 2009). The second set of Jup 



1 — 3 points 



are i nterferometr ic measurements integrated over the central 4" 
from Krips et al.l (2011). and the fainter CO(l-O) flux is the same 
from ^chmnerer et al. (2000). The shaded region shows the range 
of good-fitting LVG models (x^ < 3/2; section I3.2.3II . and the 
solid black curve is the single best fitting model, with the red and 
blue curves showing the individual contributions from the ME and 
HE components. Middle: Central velocities of the FIR CO lines, 
with average values of the ME and HE line centers indicated with 
horizontal red and blue lines. Bottom: Measured FWHM of the 
FIR CO lines, with tracks indicating the expected line widths for 
sources with intrinsic widths of 0, 100, 200, and 400 kms~^. 

In the top panel of Figure [2] we show the hne fluxes 
and upper limits measured here, along with lower-J mea- 
surements obtained from the literature. The middle and 
bottom panels show the central velocities and widths ob- 
tained from the Gaussian fitting. The inflection point 
seen in the FIR CO line SED at Juppcr ~ 19 suggests 
the presence of multiple components, as does the shift 
in central velocities between the lowest- and highest- J 
transitions. For simplicity we assume the FIR CO lines 
are produced by 2 discrete components: a moderate ex- 
citation (ME) component near the systematic velocity, 
and a blueshifted high excitation (HE) component. Our 
excitation analysis described below indicates that the 



< 17 and Juppcr > 20 transitions are dominated 



by the ME and HE components, respectively. Separately 
averaging the central velocities of these two sets of lines 
gives VyiE. = 17 ± 12 km s~^ and Vhe — — 54 ± 21 km s~^ 
(Figure [2]), with a difference of Vme — V^he — 71 ± 25 
kms^"'^. 

3.2. LVG Modeling 

To quantitatively analyze the FIR CO line SED we em- 
ploy a large velocity gradient ( LVG) model. We use the 
LVG calculation described in iHailev-Dunsheath et alJ 
(|2008D, with up d ated CO-H2 collisional coefficients 
from lYang et al] ()2010D . and a thermalized H2 or- 
tho/para ratio. In this model the shape of the CO line 
SED is determined by the gas density (nH2)i kinetic tem- 



perature (Tkin), and CO abundance per velocity gradient 
([C0/H2]/((iu/dr)). The source is assumed to consist of 
a number of unresolved clouds, and the absolute line lu- 
minosities scale with the total CO mass (Afco)- In the 
following analysis we use a CO abundance of [CO/H2] 
— 10^* to reparameterize [GO 1^2] I {dv / dr) as dv/dr, 
and Mqo as the H2 mass {Mh2)- This results in an 8 
parameter model, with 4 parameters for each of the ME 
and HE components. 

3.2.1. Background Radiation 

The CO excitation may be affected by the background 
radiation, and we must therefore estimate the local ra- 
diation density. At millimeter wavelengths the back- 
ground is domi nated by the cosmic m icrowave back- 
ground (CMB) (jKamenetzkv et al.|[20Tll) . but the FIR 
background arises from within the galaxy. The PACS 
integral field spectra provide a continuum map of the 
central 47" x 47" with 9.4" pixels, and demonstrate that 
the continuum emission in the central pixel is dominated 
by sources within the central salO". The flux density 
in the central pixel may be modeled as an optically thin 
modified blackbody with (3 = 1.5, Tdust = 48 K, and nor- 
malized to Fy{\ — 100 ^m) — 49 Jy, and we adopt this 
as an estimate of the continuum brightness of the ~5" 
CND. If the flux in the central spaxel is indeed due solely 
to the CND than this is a moderate (by a factor of < 2) 
underestimate, while emission from within the central 
«10" but outside of the CND may also be contributing. 
The measured continuum is in reasonable agreement with 
previous calculations and observations. For comparison, 
this measured Fi, is a factor of 1 — 4 times larger in 
the A = 50 — > 200 /im rang e than obtained f rom t he ra- 
diative transfer modeling in iSpinoglio et al.l (|2005D . Ad- 
ditionally, extrapolating our modeled SED to A = 450 
/im yields F^{\ — 450 /^m) = 1.2 Jy, comparable to the 
peak valu e of ^1.5 Jybeam~ ^ measu red in a ~9" beam 
by Papadopoulos fc Seaauis^ (|1999aD . The 60 /im/lOO 
/im ratio in this model is 1.27, and the FIR AujO is 
2.7 X 10-^2 Wm-2, giving Lfir = 47rd2jFiR = 1.7 x 10^° 

We include the effects of background radiation 
on the equations of statis tical equilibrium follow- 
ing |Poelma^|^iS2aani (|2005[ ). The important param- 
eter in this approach is the mean specific intensity of the 
external radiation field at the cloud surface, which we 
define as Ju,cxt- We estimate Ju,c^t using a simple ge- 
ometrical model in which the gas clouds are uniformly 
distributed in a sphere with an observed angular size 51, 
and are evenly mixed with the FIR-emitting dust grains. 
For optically thin continuum, the mean value of Ju,cxt 
is then related to the observed continuum fiux density 
-^i/,obs as 



J. 



= J 



9 F^, 



obs 



16 n 



(1) 



where /,/,cb is the sum of the CMB and cosmic IR back- 
ground (CIB), and we take 57 to correspond to a circular 
diameter of 4". We have run calculations both including 
and ignoring the local contributions to the background, 
and see negligible difference in the results. In part this 



10 Fp 



1.26 X 10-14 [2.58/6o/Jy + /loo/jy] Wm- 
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is due to the fact that the background radiation tem- 
peratures are only Trad = 13 — 21 K at the wavelengths 
of the detected FIR lines, while the typical excitation 
temperatures achieved through collisions alone for the 
best-fitting models are Tox ~ 100 K and Tox ~ 400 K for 
the ME and HE transitions, respectively. In addition, 
most of the lines are optically thick for the best-fitting 
models, and hence the CO is insulated from the external 
radiation field. 

3.2.2. Parameter Limits 

We explore two component fits to the FIR CO emis- 
sion over a large volume of 8-dimensional parameter 
space, applying physical limits to the model parame- 
ters. The most important prior restrictions are placed 
on the velocity gradient. For self-gravitating clouds 
in virial equilibrium, we can approximate {dv/dr)^i^ w 
10 kms^ipc-i (nH2/10^ cm-^f'^ (|Coldsmithl l2?f)l . 
The actual velocity gradient may be larger due to 
additional sources of gravitational potential, a high 
press ure inter-cloud inediu m, or non-virialized mo- 
tion (|Brvant fc Scovillill996D , but smaller values are un- 
likely. Defining Kyi-^ as the ratio between dv/dr and 
{dv/dr)vir ( Papado poulos et al. 2007,): 

dv/dr f nH2 V^^^ fr,. 
™ " 10 kms-ipc-i ^105 cm-3 J ' ^ ' 

we restrict parameter space to Xvir > 1- The largest 
measured velocity gradient in NGC 1068 is in the H2O 
maser disk associated with the AGN. The line of sight 
velocities of the mas er spots shift by ~600 kms~^ over a 
~2 pc linear range (jGallimore et al.ll200lD . correspond- 
ing to an effective dv/dr ^ 300 kms~^pc~^. To ac- 
commodate the maser disk and other high dispersion 
regions in our models, we extend our calculations up 
to dv/dr — 1000 kms~^pc~^. We note that restrict- 
ing i^vir > 1 and dv/dr < 1000 kms~^pc~^ combine to 
limit the density to < 10^ cm~'^, but as we discuss 
below, such high densities are ruled out by other consid- 
erations. We calculate the total gas mass in our models 
as Mgas = 1.36 X Mh-), in c ludin g the contribution from 
helium. ISchinnerer et al.l ()2000l ) estimate a dynamical 
mass of Mdyn = 9 x 10* Mq for the CND, and we dis- 
card any of our models in which the total Mgas of the 
two components exceeds Mdyn. 

The range of parameter space allowed by the CO data 
can be further reduced by considering the H2 pure rota- 
tional lines, which arise from states with similar upper 
energy levels as the FIR CO transitions. We calculate 
the H2 rotational spectrum for each model, under the 
simplifying assumption that the lines are optically thin 
and thermalized, and the H2 ortho/para ratio is thermal- 
ized. We then rule out any model that overpredicts the 
flux in any of the lin es measured in th e large (14" — 27") 
ISO-SWS apertures ()Lutz et al.ll2000[ ). These prior con- 
straints on the LVG model parameters are summarized 
in Tabled 

3.2.3. General Features of Good-Fitting Models 

We proceed by generating model SEDs over a regular 
8-dimensional grid, and for each model calculating in 



TABLE 2 
LVG Model Restrictions 

1) K^ir > 1 

2) dv/dr < 1000 kms-^pc"! 

3) 1.36 X [Mj^2(ME) + Mj^2(HE)] < 9 x 10** Mq 

4) H2 rotational lines not overproduced 



the normal manner. With 10 data points and 8 free pa- 
rameters, our modeling has 2 degrees of freedom (dof), 
and here we discuss the general properties of the set of 
solutions satisfying /doi < 3/2. In Figure [5] we show 
the range of SEDs covered by this set of good solutions, 
and for the single best flt model we show the decompo- 
sition into the ME and HE components. The 00(18-17) 
and 00(19-18) lines typically receive comparable con- 
tributions from the ME and HE components, while the 
lower- and higher- J transitions are dominated by the 
ME and HE components, respectively. The shape of the 
ME SED is relatively well constrained, and peaks in the 
>Aippor = 13—16 range. The single dish measurements 
of CO(l-O), 00(2-1), CO(3-2), and 00(4-3) we show in 
Figure [2] were obtained with 11" — 21" beams, in some 
cases comparable to the Herschel-PACS resolution (jlsraell 
12009^. However, these low- J lines receive strong contri- 
butions from the lower excitation gas in the sal5" radius 
starburst ring, and do not constrain our models. The 
interferometric measurements of the 00(1-0) flux in the 
CND range from 20 — 120 Jykms~^ (jSchinnerer et al.) 
120001 iKrips et al.l 1201 1[ ). larger than the median ME 
model flux of 8 Jykms""'^, suggesting the ME compo- 
nent contributes no more than a minor fraction of the 
observed low- J emission. The HE component is much 
less well constrained, with the best fitting SEDs that 
peak at Juppor ~ 23 only providing a marginally lower 
than those that turn over at Juppor > 30. For these 
latter models, our upper limits to 00(30-29) and 00(34- 
33) become useful constraints. 

The FIR CO emission is an important coolant of the 
nuclear molecular ISM, but does not dominate. The 
total luminosity emitted in the 10 transitions detected 
here is ico,FiR = 3.1 x 10^ Lq, and summing the mod- 
eled ME and HE emission over all transitions yields 
ico.ME + ico,HE = (4.3 - 10.5) X 10*^ Lq. For the 
highest excitation models some additional cooling may 
arise from the Juppor > 40 transitions not included in our 
LVG calculation, and significant emission is also expected 
in the Juppor < 13 submillimeter transitions. The total 
emission in the H2 0-0 S(l), S(3), S(4), S(5), and S(7) 
rota tional l ines detected by ISO-SWS is L-a^ = 1.7 x 10^ 
Lq (jLutz et al. 2000). Treating the upper limits to S(0), 
S(2), and S(9) as detections increases this by a factor of 
1.5, while at the same time some fraction of the lowest- J 
emission measured with the largest apertures may arise 
from the starburst ring. Our PACS scans have also de- 
tected a number of OH and H2O transitions. The bulk 
of the emission in these molecules detected in the cen- 
tral spaxel likely arises from the unresolved CND, and 
with this assumption we estimate nuclear luminosities of 
ioH,FiR ~ 1-5 x 10'' Lq and Lh20,fir ~ 3.0 x 10^ Lq. 
The FIR range includes the strongest OH lines at 79 ^m, 
119 /im, and 163 /im, while for H2O (as with CO) the 
longer wavelength emission should be strong. In the FIR 
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the CO and H2O cooling is therefore comparable, while 
the FIR CO luminosity is weaker by a moderate factor 
than the OH and H2 rotational emission. 

3.2.4. Bayesian Analysis 

We follow the Bayesian approach outlined 
bv lWard et al.l ()2003[ ) to quantify the probable values of 
the model parameters. We consider an 8-dimensional 
array of bins centered on the grid points for which 
we have generated a model SED, and calculated a 
value. The probability that the actual parameter set 
falls within a given bin is proportional to the product 
of bin size, likelihood L cx exp(— x^), and an assumed 
prior probability. We choose priors that are flat in the 
logarithm of each parameter, and that go to zero for 
any model that violates one of the restrictions listed in 
Table m 

In Figure [3] we show the joint density-temperature 
probability density functions for both components, with 
contours at 68%, 95%, and 99% of the enclosed proba- 
bility. The close similarity between the 95% and 99% 
(and in some regions also the 68%) contours is due to a 
truncation of the density function following the violation 
of one of our model restrictions. Starting with the ME 
component, the behaviour of the contours may be un- 
derstood as follows. For either a fixed iiTvir or dv/dr, the 
shape of the CO SED is approximately conserved if an 
increase in density is matched by an approriate decrease 
in temperature. The pair of blue curves in Figure [3] show 
the density-temperature relation best fitting the data for 
ifvir = 1 and dv/dr = 1000 kms~^pc~^. As the shape 
of the ME SED is well constrained, the region of accept- 
able parameter space for either set of models corresponds 
to a narrow band centered on the best fit curve. These 
two bands intersect at nH2 = 10^ cm~^ and diverge at 
lower densities, thereby bracketing a region of high prob- 
ability filled by models with intermediate velocity gradi- 
ents. For lower temperatures a decreasing fraction of the 
CO is excited to the higher-J states, and a larger to- 
tal mass is needed to reproduce the absolute line fluxes. 
For nH2 ^ 10^ cm^'^ the gas mass exceeds the dynamical 
mass, and consequently a small region of parameter space 
is excluded. For higher temperatures the model emission 
does not fall off with increasing J in the Jupper ~ 18 re- 
gion as rapidly as the data require, and the lower quality 
fits limit the 68% confidence interval to Tkin < 230 K. 
For Tkin > 400 and ^ 10^ cm~^ the models begin 
to overproduce the detected S(l) and upper hmit to S(2) 
H2 rotational lines, and this accounts for the sharp cutoiT 
in the probability density in the upper left region. 

The shape of the HE SED is not as well constrained 
as the ME SED, and the H2 rotational lines provide a 
stronger constraint to the extent of the probability den- 
sity contours in Figure |31 The best fit HE SEDs peak 
at Jupper ~ 23, and while lower excitation solutions are 
much less probable, higher pressure models peaking at 
Jupper > 30 also provide reasonable fits. As a result, the 
dv/dr < 1000 kms~^ pc~^ restriction no longer forms the 
upper right boundary of the density function, and mod- 
els are allowed to extend into this higher density and 
temperature region until the modeled H2 emission ex- 
ceeds the observations. The measured S(0) upper limit 
provides the most important constraint at high densities 
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Fig. 3. — Top: Joint density- temperature probability density 
function for the ME component. Contours are drawn at 68%, 95%, 
and 99% enclosed probability. The mean of log (M// 2/^^©) needed 
to reproduce the absolute line fluxes is shown by the red curves, 
and the logarithm of the thermal pressure (log(nH2 X Ty^in) in K 
cm~'^) is shown by the dotted black lines. Blue curves show the 
density-temperature relation giving the best solution for -ffvir = 1 
and dv/dr = 10^ kms~-'^pc~^ (see text). Bottom: Same as the 
top, but for the HE component. 

{nH2 ~ 10*'^ cm^'^) and low temperatures (Tkin ~ 200 
K). The higher excitation H2 lines become more impor- 
tant at higher temperatures, and the S(7) transition lim- 
its the upper left region at T^n ~ 1000 K. 

In Figure [4] we show the fully marginalized distri- 
bution functions for each of the 4 primary model pa- 
rameters, as well as for JCvir and the thermal pressure 
P/ks — X Tkin- The distribution functions for the 
HE Tkin, nH2, and P/ks are shifted to higher values 
than for the corresponding ME parameters, although 
the density distributions for the two components con- 
tain significant overlap. The lower limit imposed on K^ir 
translates into a lower limit to dv/dr that increases with 
density (see equation ^ . As such, lower density solu- 
tions are found over a broad range of velocity gradients, 
while higher density models are limited to larger values of 
dv/dr. This accounts for the positive slope of the dv/dr 
distribution. Similarly, the upper limit placed on dv/dr 
generates a negative slope in the JCvir distribution. As 
the distribution functions for these two parameters are 
heavily influenced by our prior restrictions, our model- 
ing produces no meaningful constraint on the dynamical 
state of the gas. 

We take the median of each distribution as the sin- 
gle best estimate of the parameter value, and indicate 
this number with a vertical line in each panel of Fig- 
ure m Most of the distribution functions are relatively 
symmetric, and we calculate the equivalent la uncertain- 
ties in the parameter estimation by finding the range of 
values symmetrically enclosing 68% of the total probabil- 
ity. This range is shown as the hatched area under each 
distribution. The results of this analysis are summarized 
in Table [31 The asymmetry in the ME Tkin distribution 
results from the limit to low temperature, high density 
parameter space following the K^iy > 1 and Mgas < Mdyn 



FIR CO in NGC 1068 



7 




Kinetic Temperoture [K] 




• 0.12 
0.10 
0.08 
0.06 
0.04 



>, 0.08 
I 0.06 
° 0.04 
S 0.02 
°- 0.00 




10' 10' 10" 10" 

Thermo! Pressure [K cm"'] 




10 100 

Velocity Grodient [km s"' pc"'] 




10' 10' 
Hj Moss [M,„ 

Fig. 4. — Probability density functions for the 4 primary model 
parameters, as well as for K^[^ and the thermal pressure P/kg = 
"H2 ^ ^kin- In each panel the ME distribution is shown in red, and 
the HE in blue. The binwidths are proportional to the logarithm 
of the model parameter, and the ordinate indicates the probabil- 
ity that the actual value lies in the given bin. The vertical lines 
indicate the median of each distribution, and the hatched regions 
indicate the symmetric 68% confidence interval. 

restrictions (the lower right region in the top panel of Fig- 
ure [3]). For this parameter, as well as for dv/dr and ifvir, 
we extend the acceptable range listed in Table |3] to the 
truncated edge of the distribution. 

4. COMPARISON WITH OTHER MOLECULAR TRACERS 

4.1. Physical Parameters and Mass Fraction 

The bulk of the emission traced by high reso- 
lution millimeter and near-IR molecular gas maps 
in the central 10" arises fr o m the ^5" (~350 pc 
CND (ISchinnerer et alJ I2OOOI: IGalliano fc Alloh] 12007. 
iGarci'a-Burillo et al.ll2010i r Our LVG modeling indicates 
that no more than half of the CO (1-0) emission in the 
CND is generated by our ME and HE components, indi- 
cating lower excitation material is present. The physical 
conditions of this low ex citation component h ave been 
studied by many groups. iTacconi et al.l ()1994[ ) detected 
strong CO(4-3) emission toward the center of NGC 1068, 
and combined this with an interferometric CO(l-O) mea- 
surement to show that the gas was both warm (Tkin > 70 
K) and dense {nH2 > 2 x 10"* cm~'^). Subsequent mod- 
eling of Juppor < 4 transitions of ^^CO and ^'^CO have 
typically adopted a fixed Tkin = 50 K, and derived riHn ^ 
10^ - 10^ cm-3 (iSternberg et a l.' 1991; 'Hei fer fc Bliti 
[19951 lUsero et al.l l2004D . iKrips et al.. (2011,) have ob- 



tained fluxes of the lowest 3 transitions of both ^^CO 
and ^'^CO with < 2" resolution. For a subsection of the 
CND in which the gas appears perturb ed by the radio 
jet, and hence possibly shock-heated. Kr ips et al.l ()2011l ) 
derive Tkin > 200 K and nH2 = lO^-^-'*-^ cm'^. While 
this section may not be representative of the CND as a 
whole, this analysis does suggest that the global temper- 
atures may be higher than assumed by previous authors, 
and similar to that derived for ou r ME c omponent. We 
also note that iKamenetzkv et all ()2011[ ) have similarly 
derived a globally high temperature (T > 100 K) for the 
CND based on an analysis of CS and other high den- 
sity tracers. Given this range of temperature estimates, 
the clearest difference between the ME CO component 
and the lower excitation material traced by the low- J 
CO lines is the higher density (ri//^ ~ 10^'^ cm~^ vs 
^ 10^ cm~^). This suggests a scenario in which 
the FIR lines trace denser material in the CND, which 
coexists with a more diffuse medium that generates the 
millimeter CO emission. 

The mass fraction of the high excitation gas may be 
estimated by comparing the ME and HE masses derived 
here with the mass traced by CO(l-O). With a stan- 
dard Galactic conversion factor N{H2)/Ico = 2 x 10^° 
cm~^ (Kkms~^)~^, ISchinnerer et al.l ()2000D estimate a 
mass of Mi/2 = 5 X 10^ Mq f or the CND. This is s imi- 
lar to previous estimates from iPlanesas et al.l (|199lL af- 
ter correct ing their mass to th e d = 14.4 Mpc used 
here) and iHelfer fc Bliti ()1995D . b oth of whic h used 
similar conversion factors. However, lUsero et al.l (|2004D 
have derived a much lower N{H2) / Ico = 0.3 x 10^" 
cm-2 (Kkms-i)-i (for the [CO/H2] = 10^^ used here) 
from their excitation modeling of the low- J CO emis - 
sion from t he CN D. iPapadopoulos fc Seaauisti ()1999b[ ) 
and ' Israeli (|2009[ ) have also derived conversion factors 
significantly lower than the Galactic value by analyzing 
the low-J CO emission averaged over beam sizes of ~1' 
and ~21", respectively. These authors have suggested 
the lower conversion factor arises fr om the gas n ot be- 
ing virialized. At the same time, iKrips et al.l ()2011[ ) 
have reported a CO(l-O) flux from the central 4" of 120 
Jykms"^, much larger than the 2 Jykms"^ value re- 
ported bv ISchinnerer et al.l (|20Q0l ). For N{H2)/Ico = 
0.3 X 10^0 cm-2 (Kkms-i)-i, and Fco(i-o) = 20 - 120 
Jykms~\ we estimate Mh2 = (0.5-3) x lO'^ Mq. This 
is comparable to the Mh2 = (0.1 — 8) x 10^ Mq range we 
associate with our ME component. While the uncertain- 
ties of both numbers are high, this suggests that the ME 
component makes a non-negligible contribution to the 
total mass budget. The HE emission traces a lower mass 
of warmer gas, albeit still cooler than the small amount 
(M ~ 10^ Mq) o f hot (T r^ 2000 K) gas detected in the 



near-IR H2 lines (iRotaciuc et aLl figgiTlBlietz e rallll994l : 
iGalhano fc Alloinll2002D . 

4.2. Line Profiles 

4.2.1. PACS Lines 

Further insight into the nature of the ME and HE com- 
ponents may be obtained by comparing the FIR CO line 
profiles with those of other molecular tracers. In par- 
ticular, we seek to understand the physical origins of 
the velocity shift between the two components. To bet- 
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TABLE 3 
LVG Model Results 



Parameter 


Median 


ME 

68% Range 


Median 


HE 

68% Range 




146 


7ia 


- 315 


435 


234 


- 791 


nH2 [cm ^] 


105.7 


104.9 


- lO" * 


106.3 


105.7 


- lO^-a 


P/fcs [Kcm-3] 




107-4 


- 108-« 


109.0 


108.5 


- 109-7 


dv/dr [kms~^pc~^] 


141 


26 - 


lOOO'' 


236 


56 - 


1000=* 




4.1 


la 


- 17 


3.4 


la 


- 11 


Mh2 [Mq] 


106.9 


106.1 


- lo^-'J 


105.8 


105.0 


- 106-7 



''Extended beyond the formal 68% confidence interval to the truncated edge of the distribution. 



ter demonstrate the spectral shift, in Figure [5] we sho-w 
the composite spectra obtained by averaging the ME 
{J^pp„ = 14-17) and HE (Juppcr = 20-22, 24) profiles. 
The PACS spectral resolution changes from 191 — 246 
kms~^ and 278 — 311 kms^^ over the -wavelength range 
corresponding to the ME and HE transitions, respec- 
tively. The composite spectra have been obtained by 
smoothing each line to 311 kms~^ resolution (each line 
measured "with instrumental resolution Sv is smoo thed 
with a Gaussian kernel of FWHM = v'3112 -|- Sv"^) and 
resampling to a common velocity grid. Gaussian fits 
to the composite profiles yields similar central velocities 
(^ME = 16 kms~^ and Vhe = —56 kms~^) as obtained 
from averaging the central velocities of the individual 
lines (Vme = 17 kms~^ and Vhe = —54 kms~^; Fig- 
ure [5]). The uncertainties in the centroids of the com- 
posite spectra are dominated by the spectral calibration 
uncertainties in the individual lines. With each line hav- 
ing a calibration error of ai, we estimate the error (cr) 
in the composite spectra as —< of > / N . This gives 
cme = 11 kms^^ and cthe = 14 kms~^, and we show 
these as horizontal error bars in Figure [5l 

In addition to CO, the PACS scans detect strong 
molecular emission from several transitions of OH and 
H2O, and weaker emission from 0H+, H2O+, and other 
molecules. Some of the OH and H2O lines show extended 
emission over the PACS array, but for each line the flux 
in the central spaxel is dominated by a compact source 
that we associate with the CND. The 0H+ and H2O+ 
emission is also unresolved, although an association of 
these lines with the same molecular gas is less certain. 
In Figure |6] we compare the central velocities and widths 
of these lines. For OH, each point represents the aver- 
age properties of a doublet, while for 0H+ and H2O+ 
each point is an average of multiple fine-structure tran- 
sitions. On the right hand side of the top panel we show 
the mean central velocity of each molecule, with the ME 
and HE CO components separately. The HE CO is the 
only tracer systematically offset from the systemic veloc- 
ity. In the bottom panel we show the widths of the ME 
and HE composite profiles as open points at A = 109 /im, 
where the smoothed 311 kms~^ resolution is equal to the 
instrumental resolution. Overplotted are a set of curves 
showing the expected measured linewidths for intrinsic 
linewidths of 0, 250, and 400 kms~^. The individual 
CO lines, as well as the composite profiles, are consis- 
tent with an intrinsic FWHM ~ 250 kms~^ for both the 
ME and HE components. The other molecular lines are 
somewhat narrower. 

The ^50 — 70 kms^^ offset between the HE lines and 
the other molecular tracers in Figure [6] is <l/3 of the 
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Fig. 5. — Average profiles of the ME (red) and HE (blue) CO 
lines. Prior to averaging, each line is smoothed to a resolution 
of 311 kms~^. The horizontal error bars shown underneath the 
centroid of each profile indicate the estimated ilcr calibration un- 
certainty of the stacked spectra. 

PACS spectral resolution, and at this level instrumental 
effects must be considered. The PACS spectral response 
depends on the illumination of the slit, such that a phys- 
ical translation of a source on the sky in the dispersion 
directio n will produce a wave length shift in the spectral 
profile (jPoglitsch et al.|[2010D . In the top panel of Fig- 
ure [6] we show the equivalent velocity shift resulting from 
moving a point source ±2" from the center of the slit. A 
comparison with the CO velocities shows that a ~2 — 3" 
offset between the centers of the slit and the HE CO 
emission could be partially responsible for the measured 
wavelength shifts. With the exception of CO(17-16) and 
00(24-23), each point in Figure [6] corresponds to a line 
measured in the same set of concatenated range scans. 
Therefore if all of the molecidar emission is presumed to 
arise from the same region on the sky, a wavelength shift 
in the HE CO lines induced by a pointing offset should 
be matched by a comparable shift in the A « 100 — 150 
/zm OH and H2O lines, and this is not observed. Al- 
ternatively, the HE CO lines may be modeled as arising 
from a region centered ^2 — 3" away from the source of 
the OH and H2O emission. However, the 00(24-23) line 
at A = 109 /im was observed 6 months prior to the se- 
ries of concatenated scans, and the slit position angles of 
these two epochs differ by «180°. Any pointing- induced 
shift to the 00(24-23) line should then have a compa- 
rable magnitude but be in the opposite direction of the 
shifts in the other HE CO lines, and this is also not the 
case. We therefore argue that the observed wavelength 
shift represents a real velocity offset, with instrumental 
effects playing no more than a minor role. 
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Fig. 6. — Top: Measured line centroids for the CO, OH (green), 
H2O (brown), OH+ (purple), and H2O+ (cyan) lines detected in 
our PACS range scans. The CO lines are color coded by ME (blue), 
HE (red), and intermediate (black). Each OH point is an average 
of a doublet, and the OH"'' and H2O"'" points represent averages 
over all detected fine-structure lines. The dashed line shows the 
instrumental wavelength shift induced by pointing offsets of ±2" in 
the dispersion direction. On the right hand side we show the mean 
centroid of each tracer. Bottom: Same as the top, but showing 
the measured FWHM of each line. At A = 109 fim we show the 
linewidths of the ME and HE composite profiles as open symbols. 
Solid curves show the expected measured widths for lines with 
intrinsic widths of 0, 250, and 400 kms~^. 

4.2.2. High Resolution Observations 

To search for kinematic substructure in the CND that 
might explain the velocity shifts among the FIR CO lines, 
we compare the FIR line profiles with those of high res- 
olution observatio ns of C0(2-l) (ISchinnerer et al.N2000D 
and H2 1-0 S(l) (|Miiller Sanchez et all 120091) . In Fig- 
ure [7^ we show maps of these two tracers in the central 
4" X 6". At 0.7" resolution the C0(2-l) emission from 
the CND is dominated by a pair of knots centered ~1" 
east and r^l.5" west of the AGN, connected by lower 
surface brightness emission. Similar morphology is also 
seen in th e ^1" resolution aperture sy nthesis maps of CN 
and Si O (iGarcia-Burillo et afl [2010I ) and of HON and 
HCO+ (|Krips et all 1201 It) . Observations of the H2 1- 
S(l) line at 0.075" resolution have resolved the CND 
into a r^l.2" radius ring (hereafter the "H2 ring") cen- 
tered ~0.6" southwest of the AGN. The eastern knot is 
also prominent in the II2 map, while the western knot is 
much fainter. The H2 map also shows strong emission 
from a clump ^1" to the north of the AGN that is not 
prominent in the C0(2-l) map. 



In panel b we show separate maps of the C0(2-l) 
emission int egrated over blue an d red velocities (see also 
Figure 6 in iKrips et all ()2011l) ). The strong emission 
to the east is largely to the blue of the systemic ve- 
locity, while the emission to the west separates into a 
blueshifted SW component and a redshifted NW com- 
ponent. The H2 velocity field is generally similar, with 
the exception of t he no rthern clump as discussed below. 
ISchinnerer et all (|2000f ) model the C0(2-l) kinematics 
as a warped disk, while subsequent study of the H2 and 
millimeter tracers has shown evidence for additional non- 
circular mot ion, possibly indicating a radially expanding 
component (iGalliano fc AlloinI 120021: iDavieset al.l 120081: 
IGarcia-Burillo et all 120101: iKrips etal.ll2011[ ). hTpanels 
e through i we compare the composite ME and HE spec- 
tral profiles with those of C0(2-l) and II2 extracted from 
selected apertures. For each panel we have smoothed 
the C0(2-l) and II2 spectra to the same resolution (311 
kms~^) and resampled to the same grid as used for the 
ME and HE composite spectra. 

The most important result of this comparison is the 
mismatch between the FIR line profiles and the broad 
and highly blueshifted H2 emission from the bright 
knot ^l" to the north of the AGN, shown in panel 
f. In the high spec t ral re solution H2 map presented 
in IGalliano fc AlloinI ()2002t ) this is the one region in 
the H2 ring that displays a double-peaked profile (the 
two peaks become blended following the smoothing done 
here), with an extra emission component to the blue that 
is a clear outlier to their simple rotation plus expansion 
model. This region is also the site of the strongest Br7 
emission in the H2 ring, which displays a broad (FWHM 
^ 1000 k ms~^) line. The radio jet exits the nucleus to 
the north (iGallimore et al.lll996( ). an d may interact wit h 
material in the narrow -l ine region (lAxon etlII[T998h . 
IGalliano fc AlloinI (|2002| ): lGamano et al.l(|2003f ) interpret 
the broad and complex Br7 and H2 line profiles as arising 
from gas perturbed by the jet, and additionally suggest 
that the ionization resulting from the jet-ISM interaction 
may be responsible for the strength of the Br7 emission in 
this region. The smoothed H2 profile in panel f is incon- 
sistent with either the ME or HE composite spectrum, 
and we conclude that this region of the CND does not 
dominate the PACS CO emission. In sections |S] and [51 
we discuss possible excitation mechanisms for the FIR 
CO, including shock heating. The fact that we can rule 
out an origin in the region of the H2 ring with the best 
evidence for molecular gas perturbed by the jet leads us 
to suggest that such jet-driven shocks in the H2 ring do 
not excite the high- J CO. 

Aside from the bright H2 knot in the north, the ME 
profile is generally consistent with much of the rest of the 
H2 ring. The best fit is with the profile of the strong H2 
emission from the east. The profiles from the NW and 
SW are moderately red and blueshifted with respect to 
the ME composite, but a combination of these regions, 
and indeed of the H2 emission integrated over the entire 
CND (excluding the northern knot), would also gener- 
ate a reasonable fit. We note that the H2 linewidths 
(FWHM = 390 - 440 kms"!) are larger than the C0(2- 
1) linewidths (FWHM = 330-350 kms'^) in each panel, 
and better match the widths of the composite ME and 
HE profiles (FWHM = 400 - 410 kms^i). This may 
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Fig. 7. — a) CO(2-l) (contours) and H2 1-0 S(l) (color) images of the central 4" X 6" from lSchinnerer et"an l|2000l ): IMtiller Sanch ez et al] 
1(2009). b) blue and red components of CO(2-l), with crosses marking the centers of the apertures in a), c) H2 1-0 S(l) image of the central 
0.4", e-i) smoothed H2 1-0 S(l) (black solid) and CO(2-l) (black dashed) profiles from selected apertures in the H2 ring compared with 
the ME (red) and HE (blue) composite profiles, and ±40 kms~^ horizontal error bar indicating the calibration uncertainty between the 
FIR CO and H2 spectra, d,j) smoothed H2 1-0 S(l) profiles from the northern and southern streamer compared with the ME and HE 
composite profiles, and ±40 kms^^ horizontal error bar. 



indicate that the hot gas probed by the H2 is a better 
tracer of the material producing the high- J CO, although 
due t o extinction of the H2 hne (see, e.g., [Galliano et al.l 
l2003f ). the morphology of the FIR CO emission may be 
different than the H2 image. 

The blueshift of the HE emission is more challenging 
to match. The H2 spectra from the blue regions in the E 
and SW are «60 and «4G kms"^ to the red, respectively, 
while the CO (2-1) profiles from the same regions are too 
narrow. However, the uncertainty in the mean centroid 



of the HE lines is «20 kms^^, while the H2 calibration 
error is smaller. Registration errors between the veloc- 
ity frames of the PACS spectra and the ground-based 
H2 spectra are also likely to be important at the level 
of ~10 — 20 kms~^, although more difficult to quantify. 
With a conservative estimate of 40 kms^^ uncertainty 
in the relative calibration of the PACS CO and the H2 
spectra, the HE CO emission only differs by 1 — 1.5a 
with the H2 emission from the E or SW regions. We 
conclude that while the HE emission profile is not natu- 
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rally matched to the H2 or C0(2-l) emission from the H2 
ring, an association with the bluest material to the E or 
SW may be within the measurement errors, and should 
not be excluded. 

In panel c we show the 0.025" resolution H2 1-0 S(l) 
map, which identifies two gas clouds streaming toward 
the nucleus on highly elliptical orbits from the north 
and south (hereafter the "northern" and "southern" 
streamers; iMiiller Sanchez et al.l I2OO90 . The northern 
streamer is connected to the II2 ring in both H2 emis- 
sion and mid-IR continuum, and has been proposed 
as a r neans by which material is transported to th e 
AGN (|Tomono et al.l[200l iMOller Sanchez eTall 120091) . 
The southern streamer is detected to within ~10 pc of 
the AGN, and may play a role in obscuring the nuclear 
emission. The southern streamer is modeled to lie in 
front of the AGN in the plane of the galaxy, and has 
a redshifted velocity that increases from ~50 kms~^ to 
^80 kms~^ as the gas moves to within a projected dis- 
tance of < 0.1". The northern streamer approaches the 
AGN from behind, and the brightest emission occurs 0.4" 
from the center at a blueshifted velocity of ^ — 30 km s^ ^ . 
In panels d and j we compare the smoothed profiles of 
these two regions with the PACS lines. The emission 
from the southern streamer is detected with lower S/N, 
but displays a profile reasonably consistent with that of 
the ME CO composite. The line profile of the northern 
streamer produces an excellent match to the HE compos- 
ite. The centroids of the southern streamer and the HE 
CO composite differ by salOO kms~^, which we argue is 
too large to be plausibly explained by calibration uncer- 
tainties between the two datasets. We conclude that in 
addition to an origin in the H2 ring as discussed above, an 
origin of the ME and HE components with the southern 
and northern streamers, respectively, would be consistent 
with the line profiles. 

5. HEATING THE ME COMPONENT 

The ME CO emission arises from a warm (Tkin ~ 146 
K) and dense {tih^ ^ 10^-'' cm~^) component, and with a 
total mass of Mjj^ ~ 10^'^ Mq, represents an important 
fraction of the ISM in the CND. The kinematic analysis 
presented in section [4.21 shows that this component may 
be readily attributed to the H2 ring, or possibly to the 
southern streamer. Here we consider potential heating 
mechanisms, and conclude that X-ray and shock heating 
are both plausible, while heating by far-UV photons is 
less likely. We further argue that no plausible heating 
mechanism is consistent with an origin in the southern 
streamer, and hence the emission is likely associated with 
the H2 ring. 

5.1. X-Ray Heating 

The AGN in NGC 1068 emits a hard X-ray luminos- 
ity of L2-iokeV = 10"^^^"^^ ergs-i (jlwasawa et al.l 119971: 
iColbert et alll2002[ ). Our view of the AGN is obscured 
by a Compton-thick medium, but the extended emission 
detected by Chandra in the 6 — 8 keV band demon- 
strates that the n uclear X-rays irradiate the ISM over 
the c entral ~kpc ()Ogle et al.H2003t iGarcia-Burillo et "all 
I2OIOD . Hard X-rays penetrate deeply into clouds, and 
efficiently heat large columns of molecular gas through 
photoionization heating (jMalonev et al.. .1996i ). The 



bright H2 1-0 S(l) emission in the CND has been 
attributed to X-ray heated gas (iRotaciuc et al.l 119911 : 
lMalonevlll997l: iGalliano fc Alloinll2002D . and X-ray heat- 
ing s hould also produce st rong emission in the FIR CO 
lines (Krolik & Lcdd 1989). 

Mei jerink fc Spaans (2005) and lMeijerink et al.l ()2007l ) 
present a detailed photochemical modeling of X-ray dom- 
inated regions (XDRs) that includes predictions for the 
emergent CO line intensities as a function of the gas den- 
sity (tt-h) and incident hard X-ray fiux (Fx = ^2-iokcv)- 
Wc use their type A models, which calculate the emis- 
sion from a parsec thick cloud over a grid covering 
riH = 10"* -lO'^-^ cm-3 and Fx = 1.6-160 ergcm'^s'^ 
These models generate CO line SEDs with similar shapes 
as the isothermal models used in our LVG analysis, and 
an analogous two component fit reproduces the FIR CO 
line fluxes. In Figure [8^ we show a model that uses 
riH = lO'^'^^ cm"'^ and Fx = 9 ergcm~^s~^ for the 
ME emission (see Table HI . For an AGN luminosity of 
-^2-iokoV = 10^'^"'*'* ergs~^, geometric dilution of the ra- 
diation field at the d ~ 100 pc distance of the H2 ring 
yields a fiux of Fx = 8.4 — 84 crgcm^^s"^, broadly 
consistent with this mo deled flux. For the pla ne-parallel 
geometry employed by iMeiierink et al.l (|2007f ) the abso- 
lute line luminosities scale with the total XDR surface 
area, and the normalization of the ME component of 
the model shown in Figure requires A ^ (120 pc)^. 
Gallia no et al.l (|2003D model the H2 ring as a section of a 
40 pc thick disk. For a radius of 100 pc the inner surface 
area exposed to the AGN is then ^(160 pc)^, similar to 
the XDR model requirement. In sum, we conclude that if 
a substantial fraction of the H2 ring is exposed to nuclear 
hard X-rays, then both the shape of the ME segment of 
the CO SED and the absolute line fluxes are naturally 
reproduced in this XDR model. 

5.2. Far-UV Heating 

Far-UV (FUV; 6 eY < hv < 13.6 eV) photons offer 
another means of heating the molecular gas. Photodis- 
sociation regions (PDRs) powered by the FUV radiation 
from OB stars in the Galaxy have indeed been iden- 
tified a s prominent sources of FIR CO emission (e.g., 
iKramer et al.. 2004 ). as have the FUV-irradiated cavities 
of protostellar outflows (jvan Kempen et aLll2010t ). How- 
ever, a comparison with the high- J CO emission from the 
prototypical starburst galaxy M82 suggests that the CO 
emission from NGC 1068 is too strong to be powered by 
stellar FUV. Herschel observations of M82 have detected 
the submillimeter CO transitions ( Jupper = 4—13), which 
trace emission from gas h eated by shocks (iPanuzzo et al.l 
120101 ) and/or in PD Rs (|Loenen et all 12010^ . For an 
Lftr ~ 3 X 101° Lq ()Telesco fc H"a7^[T98a iJov et al.l 
Il987t ). the CO(13-12)/FIR ratio integrated over the cen- 
tral 43.4" (-821 pc) is -6 x 10"^. The FIR continuum in 
M82 is produced by the FUV output of young stars, and 
we interpret this CO(13-12)/FIR ratio as a benchmark 
estimate of the fraction of FUV radiation that may be 
converted to FIR CO emission in an extragalactic star- 
burst. The CO(14-13)/FIR ratio we estimate for NGC 
1068 is ~3 X 10-^ a factor of ~5 larger than the C0(13- 
12)/FIR ratio in M82. This discrepancy is further in- 
creased by noting that only a minor fraction of the FIR 
continuum from the CND in NGC 1068 originates in 
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TABLE 4 
Heating Mechanisms 







ME 


HE Pull 


XDR 


nn = 




riH = 10^'"' cm^'^ 




Fx = 


9 erg cm~^ 


Fx = 160 ergcm~^s~^ 




A 


~ (120 pc)2 


A~(25pc)2 


PDR 






n = 10"^-^ cm"^ n = 10" cm'-' 
Go = 10'i''5 Go = 10^ 
Lpuv ~ 2 X 109 Lq Lpuv ~ 10^" 


shock 




C-shock 


C-shock 




no = 


2 X 10^ cm-3 


no = 10^ cm~^ 




V - 


= 20 kms-i 


V = 30 km 




A 


~ (130 pc)2 


A~(21pc)2 



N ote. — Details for the mod e ls use d in Figure [S] XDR and PDR models arc from 'Meiierin k et al.l l|2007l ). ME C-shock model is 
from lFlower Ik Pineau Des Foretj 1 12013 ). and HE C-shock model is from|Kaufman & Ncufcld (199f^^ 
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Fig. 8. — FIR CO line fluxes and upper limits measured here, 
along with lower- J lines from the literature (see Figure [2] for de- 
tails), and XDR, PDR, and shock model fits, a) Overplotted is 
a two component XDR fit (solid black), with individual compo- 
nents in red and blue. Dotted section of the red curve shows an 
extr apolation of the model. Upper limits to the IKrolik He Lepd 
([1989) torus model are shown by the black (constrained by the 
CO [44-43] limit) and green (constrained by the stacked limit to 
selected Juppor = 34 — 47 transitions) dashed curves, b) Black 
curve shows a PDR fit to the full CO SED, and blue curve shows 
a separate fit to the Juppcr > 19 lines, c) Two component shock 
fit with same color scheme as panel a). Blue curve shows a model 
interpolated and extrapolated from a finite number of transitions 
(diamonds). Model parameters and references are discussed in the 
text and summarized in Table |4] 



young stars. The nuclear stellar cluster in NGC 1068 
has a characteristic age of 200 — 300 Myr, and hence 
Lbo\/LK ^ 70 (for a starburst with exponential decay 
tiniescale of tsf = 100 Myr; iD avies et al. 2007). For the 
Lk measured by IDavies et al.[ (p007[ ) this corresponds to 
-^boi ^ 3 X 10^ Lq, more than 5 times lower than the LpiR 
measured here. PDRs heated by the output of this clus- 
ter would therefore have to be ^25 times more efficient in 
converting the incident FUV radiation to FIR CO e mis- 
sion th an the PDRs in the nucleus of M82. iLoenen et al.l 
(|2010D model the Juppcr = 1 - 13 CO emission from M82 
with a 3 component PDR model, in which the C0(13- 

12) emission is almost entirely generated by the high- 
est excitation component. The CO(13-12)/FIR ratio of 
this component is comparable to the measured C0(14- 

13) /FIR ratio in NGC 1068, but the total modeled ratio 
is a factor of ~14 lower due to the FIR emission from the 
lower excitation material. Increasing the net ratio by a 
factor of ^^25 would require the total FIR be dominated 
by the highest excitation component, and hence a signif- 
icant suppression of the lower excitation PDR emission. 
We consider that such an extreme repartitioning of the 
stellar FUV between M82 and NGC 1068 less likely than 
the XDR scenario discussed in section ISTTl 

The CO emission is more plausibly attributed to PDRs 
powered by the AGN FUV radiation, and in many ways 
this is an attractive option. The intrinsic AG N luminos- 
ity in NGC 1068 of Lpuv ~ 7.4 x ergs"! (IPier et al.l 
|1994|) is similar to the observed Lfir, supporting the 
idea that the FIR continuum is emission from dust 
grains heated by the AGN. The measured 60 ^m/100 
fim flux ratio of ^1.3 implies that the incident FUV 
flux on these grains corresponds to Gp ^ 10^ (where 
i^FUV = Goxl.6xl0-3ergs-icm-2) (|Abel et al.l feoOQl. 
which would also be expected if the AGN Lpuv is ab- 
sorbed at the d ^ 100 pc distance of the H2 ring. To es- 
timate the CO emission produced in these PDRs we em- 
ploy the type A PDR models of Meijerink et al. (2007), 
which adopt the same geometry and gas densities as 
the XDR models, and span a FUV intensity range of 
Go = 10^ - 10^. In dense PDRs most of the CO forms at 
cloud depths greater than Ay « 2 at gas temperatures 
T < 100 K, but a secondary CO abundance peak at 
Av ~ 0.6 with T ^ 800 K follows from a local enhance- 
ment of OH (Sternberg fc Dalgarnolll995D. T he CO line 
SEDs generated in the lMeiierink et al.l ( 20071 ) models are 
characterized by two distinct peaks, which we associate 
with these warm and hot CO phases. In Figure [5}d we 
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show a model with nn = 10^ cm^^ and Gq — 10^ that 
produces an intriguing match to the entire FIR CO hne 
SED (sohd black line), with the warm and hot phases 
reproducing the ME and HE components, respectively. 
Assuming the emitted FIR continuum is equal to the in- 
cident FUV, this model predicts a CO(14-13)/FIR ratio 
of 4 X 10~^. This is similar to the observ ed ratio in high- 
Go Galactic PDRs ([Kramer et al.ll2004 ). and the value 
measured here for NGC 1068. 

However, there are at least 2 reasons to reject this sce- 
nario. First, the AGN UV/optical emission escapes to 
d ~ 1 00 pc distances only i n the ionization cone (PA 
w 15° IMacchetto et al.lll994D . while the brightest emis- 
sion from the H2 ring is to the east. If the warm molec- 
ular gas traced by the ME CO emission is assumed to 
be cospatial with the hot gas traced by the H2 1-0 S(l) 
line, then we would expect little interaction between this 
warm gas and the nuclear FUV. Additionally, any single 
PDR model that simultaneously matches the ME and 
HE emission is inconsistent with the kinematic evidence 
that these two sets of lines are tracing physically distinct 
components. In general, we require that any model re- 
producing the Juppcr ^17 lines must underpredict the 
fluxes in the Juppcr 20 transitions. In contrast, all of 
the models produced by iMeijerink et al.l (|2007f ) that fit 
the ME section of the CO SED either match or exceed the 
observed HE lines fluxes. Consequently we exclude PDRs 
as a potential source of the ME emission. We stress that 
this latter argument depends sensitively on an accurate 
modeling of the Juppcr > 20 CO emission from the hot 
surfaces of PDRs. Future Herschel-PACS studies of the 
FIR CO line SEDs in Galactic PDR templates will be 
useful in further evaluating this PDR scenario. 

5.3. Shock Heating 

Shock heating offers a simple means of exciting the 
high- J CO emission, and is generally the preferred mech- 
anism for producing the highest excitation molecular gas 
in Galactic sources (e.g., Scmoere et al. 2000). Using the 
C-shock models of Flower & Pi neau Des ForetsI (|2010( ). 
we find that the ME component of the CO SED can be 
fit with a preshock density of nn = 2 x 10^ cm~'^ and 
shock velocities of = 10 — 20 kms~^. The contribu- 
tions of H2O and H2 to the total cooling budget increase 
rapidly with shock velocity in these models, and keeping 
Us < 20 kms^^ is necessary to prevent the predicted H2O 
and H2 line fluxes from exceeding the measured values. 
Decreasing the shock velocity from Wg = 20 — )■ 10 kms~^ 
generates weaker CO lines, and requires increasing the 
total cross-section from A ^ (130 pc)^ — (200 pc)^ to 
match the absolute line fluxes. In Figure [5J: we combine 
the Vs = 20 kms~^ model with a separate C-shock fit to 
the HE CO emission. 

What shock mechanism could produce the ME CO 
line emission? Mechanical heating from stellar feed- 
back has been proposed as the energy source be- 
hind t he submillim eter CO emission in the M82 star- 
burst (|Pa^^z^^FaL_^010) and i n other galactic nu - 
clei (jHailev-Dunsheath et all 120081: INikola et all 1201 ID . 
but the contrast in CO/FIR ratios between NGC 1068 
and M82 discussed in section 15.21 argues against a simi- 
lar mechanism here. Furthermore, we can use the results 
of .Davies et al.. (,2007i ) to estimate the mechanical power 



injected into the ISM by supernovae (SN) in the nuclear 
cluster of NGC 1068. These authors model the super- 
novae rate (SNR) to Lk ratio as < 6 x 10~^^ yr~^ Lq^ 
in this cluster, which combined wit h their measured Lk 
yields SNR < 2.4 x 10"^. Following Eoenen et all (|2008D 
and estimating a mechanical energy release of 10^^ erg 
per SN, with 10% dissipated in molecular gas, yields a 
heating rate of <2 x 10^ Lq. This is at least a factor 
of ^8 — 28 times lower than the mechanical luminosity 
L = l/2pv^A required by the shock models discussed 
above. Jet-ISM interactions are another potential source 
of shock heating in Seyfert nuclei, but as discussed in sec- 
tion 14.21 the mismatch in line profiles between the FIR 
CO lines and the disturbed H2 1-0 S(l) profile --1" north 
of the AGN suggests jet-driven shocks in the H2 ring may 
not be important. 

Alternatively, we note that the cross-sections required 
to normalize these plane-parallel shock models are sim- 
ilar to the estimated cross section of the H2 ring as 
viewed from the center. The kinematics drawn from the 
high resolution H2 1-0 S(l) and millimeter- wave molec- 
ular gas ma ps require a radial expansion component to 
the H2 ring (iGalliano fc AUoinl [20021: iDavies et al.| [2008l: 
IKrips et al.l[2011h . IGalhano fc AUoinl (I2002D model the 
H2 dynamics by combining a rotational component with 
a V = 140 kms^^ radially expanding com ponent that 
generates 1/3 of the total line emission, and iKrips et al.l 
( 2011 ) construct a similar model to explain the CO dy- 
namics. We suggest that the interaction of the outfiow- 
ing gas with non-outflowing material in the H2 ring of- 
fers the most plausible source of shock heating. We re- 
call that our excitation model requires the dense gas as- 
sociated with the ME CO emission to be mixed with 
lower density {nH2 ^ 10^ cm~'^) material responsible for 
the millimeter- wave CO emission (section 14. ip . Assum- 
ing the nv'^ product is conserved for shocks propagatin g 
through an inhomogeneous medium (jKlein et al.lll994) . 
the moderate {v < 20 kms~^) velocities we require could 
ultimately be produced in v 140 kms~^ shocks trig- 
gered in lower density gas. 

5.4. Southern Streamer 

In section we showed that the ME CO line profile 
may be consistent with that of the H2 1-0 S(l) emission 
from the southern streamer. However, such an associa- 
tion is inconsistent with either the X-ray or shock heat- 
ing scenarios outlined above. In both of these models 
the required cross-section oi A ^ (120 pc)^ — (130 pc)^ is 
significantly larger than the ^20 pc size of the southern 
streamer. Equivalently, the absolute intensity of the CO 
emission from this cloud would have to be more than an 
order of magnitude larger than predicted by the models. 
As such, we argue that the ME CO emission does not 
arise from the southern streamer, but most likely arises 
from the H2 ring. 

5.5. Summary and Discussion 

In summary, we conclude that the ME CO emis- 
sion arises from either X-ray or shock-heated gas 
in the H2 ring. The challenge of unambiguously 
determining the heat source for this gas is sim- 
ilar to the situation for the warm and hot gas 
traced by the H2 rotational and ro-vibrational lines 
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in NGC 1068 (|Rotaciuc et all [TOQlt iLutz et alj [2OOOI) 
and in larger samples of Seyferts (iDavies et alj 120051 : 
iRodn'guez-Ardila et all [20051: iRoussel et all 120071 ) . The 
link between the ME CO and the H2 emission in NGC 
1068 is not immediately clear, although at least to within 
the PACS resolution the similarity of the ME CO and H2 
1-0 S(l) line profiles suggests a connection (Figure [7]). 
Additionally, we note that the H2 1-0 S(l) brightness 
is quantitatively reproduced with a similar XDR model 
as used here for the ME CO (although using a lower 
densi ty of n = 10^ cm~^; iMalone v 199.7,; iGalliano et alJ 
l2003f ). while the shock model shown in Figure [8]: also 
produces an H2 1-0 S(l) flux within a factor of ~2 of the 
total CND emission. Further joint modeling of the CO, 
H2, and other molecular emission in NGC 1068, and FIR 
CO data on a larger sample of comparison sources, will 
be useful in better understanding the nature of the ME 
CO component. 

Here we offer two reasons for preferring the X-ray heat- 
ing scenario. First, the hard X-ray luminosity of NGC 
1068 is reasonably well established through either an 
analysis of the directly observed (scattered) emission, or 
by ap plying scaling r elations established for ty pe 1 sys- 
tems (|Iwasawa et all 119971: IColbert et al.ll2002l) . As dis- 
cussed above, combining this luminosity with reasonable 
estimates of the gas density and the CND geometry nat- 
urally produces the ME CO line SED. In contrast, it is 
not clear whether the shocks in the CND are sufficient to 
dissipate enough mechanical energy at the low velocities 
needed to power the CO. As we argued above, jet-driven 
shocks in the II2 ring do not provide a good fit, while the 
energetics of the possible shocks arising from the radial 
expansion of the H2 ring have yet to be demonstrated. 
Secondly, a number of authors have noted that the chem- 
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[2OIO; Rangwala ct al. 2011), and will be discussed in a 
future paper. If the nuclear X-rays are responsible for 
the anomalous molecular abundances in the CND, it is 
likely they also play an important role in the energetics. 

6. HEATING THE HE COMPONENT 

The HE CO emission arises from a small mass {Mh2 ^ 
lO^-S Mq) of warm (Tki„ ~ 435 K) and dense (n//, - 
10^ '^ cm~^) material that represents only a minor frac- 
tion (< 10%) of the total gas in the CND. The kine- 
matic analysis presented in section suggests that this 
component may potentially be associated with the most 
blueshifted emission in the east or west of the H2 ring, 
or with the clump of infalling gas ^^0.4" north of the 
AGN. Here we consider potential heating mechanisms, 
and argue that the HE CO emission arises from either X- 
ray-heated gas in the northern streamer, or shock-heated 
material in either the northern streamer or H2 ring. 

6.1. X-ray Heating 

In section EU we used the lMeiierink et al.l (|2007D XDR 
models to generate a two component fit to the FIR CO 
emission, and overploted a sample solution in Figure [8^. 
For the nn = 10^ -)■ 10^-^ cm"^ density range the ME 



component requires Fx — 5.1 — )■ 16 ergcm~^s~^, while 
for the same densities the HE component requires the 
ma ximum (or highe r ) Fx = 160 ergcm~^s~^ used in 
the IMeiierink et al.l (|2007f ) model grid. An important 
outcome of this XDR modeling is therefore that the HE 
component requires irradiation by an X-ray field at least 
an order of magnitude stronger than the ME component. 
This is difficult to achieve if both components are pre- 
sumed to arise from the H2 ring. The eastern segment 
of the ring lies no more than a factor of ~1.5 closer to 
the AGN than the western segment, so a variation in the 
distance to the AGN across the ring appears unlikely to 
produce a factor of ~10 varia t ion in the radiation field 
strength. iGarcfa-Burillo et "all ()2010D have used the 6 — 8 
keV emission detected with Chandra to trace the pene- 
tration of nuclear hard X-rays into the cold neutral ISM. 
Their image suggests a relatively isotropic il lumination 
of th e central few hundred parsecs (see also lOgle et ahl 
120031) . Specifically, while the hard X -rays may be ex- 
pected to escape to ^100 pc scales more easily in the 
ionization cone, the 6 — 8 keV band image is not signifi- 
cantly enhanced along this direction within the H2 ring. 

The most straightforward way of generating a higher 
X-ray intensity is to attribute the HE emission to the 
northern streamer. In projection, the northern streamer 
lies a factor of ^3 clos er to the AGN than does t he H2 
ring (-0.4" vs ~1.2": iMimer Sanchez eralll2009l ). and 
is therefore expected to see a ~9 times stronger Fx- A 
fit to the HE component with nn = 10^'^^ cm~^ and 
Fx ~ 160 ergcm^^s^^ (Figure [SJi) requires a surface 
area of ~(25 pc)^, only a moderate factor larger than 
the ~14 pc (~0.2") projected lateral size of this clump. 
We therefore argue that if X-ray heating accounts for 
both the ME and HE components, the HE component is 
most plausibly associated with this infalling gas. 

6.2. Far- UV Heating 

Using the IMeiierink et al.l (|2007D PDR models we find 
that for a proper choice of high density (nn > 10^ cm~^) 
and strong FUV field (Go > 10^), the models reproduce 
the Juppor > 20 lines while underpredicting the lower- 
J fluxes. An example model with tih = 10^'^ cm~'^ 
and Go — lO'* '''^ is shown in Figure |H)d (thin blue line). 
The FUV continuum associated with this component is 
LpuY ~ 2 X 10^ Lq, a factor of ~10 less than the ob- 
served LpiR- Thus we may consider a scenario in which 
the ME CO emission arises from X-ray or shock-heated 
gas in the H2 ring, while the HE emission originates in a 
trace amount of PDR material that makes only a minor 
contribution to the total continuum emission. 

Any gas in the CND exposed to the FUV from the 
AGN must also be irradiated by hard X-rays. For a 
FUV origin of the HE CO lines the physical proper- 
ties and radiative environment of the clouds must there- 
fore be such that the FUV is more effective than the 
X-rays in generating Juppcr > 20 CO emission. Un- 
der what conditions would this occur? A first approach 
to addressing this question is to separately consider the 
CO emission from a PDR with that of a properly se- 
lected XDR. In Figuredlwe compare the 00(19-18) (the 
highest- J tra nsition calculated over the full model grid) 
flux from the Meijerink et al.' (2007') PDR models with 
that of a ,Meiierink ct al.. (,2007. ) model XDR with wl/6 
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Fig. 9.— Ratio of CO(19-18) emiss ion from a PD R to an XDR 
(black), using the models of Mciierin k et al.l 120071) . The incident 
X-ray flux in the XDR model is 6 times weaker than the FUV flux 
in the PDR model. The fraction of the measured LpiR that must 
be attributed to this PDR is shown in red. The region of high 
density and high-Go that can match the HE CO line SED without 
overpredicting the lower- J fluxes is indicated with the dotted lines. 



of the incident flux. This flux ratio is chosen to match 
the intrinsic L-puy/Lx ~ 6 ratio of the AGN in NGC 
1068 (|Pier et al.l 119941 ) . and is therefore appropriate for 
a cloud seeing the unattenuated AGN emission. The 
red curves in Figure [H] show the fraction of the observed 
FIR continuum (assuming Lpuv = ^fir) that must arise 
from the model PDR in order to account for the absolute 
CO(20-19) line flux. The region in the upper right quad- 
rant enclosed by the dotted lines indicates the subset of 
high excitation PDR models that reproduce the HE CO 
lines while underpredicting the lower- J fluxes - a require- 
ment for this solution given the velocity shifts between 
the ME and HE hues. For clouds with nn > 10^ cm~^, 
Go > 10^, and Go/nn < 10" cm^, this approach sug- 
gests that the FUV is more important than the X-rays in 
generating Juppcr <^ 19 CO emission, while only a minor 
fraction of the AGN FUV luminosity would be required 
to power such a PDR. 

Where in the nuclear region might these conditions 
be satisfied? The UV/optical emission from the AGN 
escapes to large distances within the ionization cone, 
which runs roughly north-south in projection. If the HE 
CO emission arises from the H2 ring, however, the line 
blueshifts would suggest an origin to the east or west 
(section 14.21) . where little nuclear FUV penetrates. The 
northern streamer may have a direct view of the AGN, 
which at a distance of d « 40 pc would correspond to 
Go ~ lO^ ''. However, the covering factor of this mate- 
rial (for a size of ~14 pc) is only ^0.01. The contents 
of Figure ini indicate that achieving the absolute HE CO 
line fluxes with such a small covering factor would re- 
quire densities far larger than the nH2_2L ^^^'^ cm^^ 
indicated by our LVG modeling (section [3. 2p . We there- 
fore argue that the HE CO emission is unlikely to be 
FUV-powered, although more detailed modeling of the 
conversion of FUV radiation to high- J CO at high den- 
sities (nn > lO''-^ cm~3) and FUV intensities (Go > 10^), 
and in the presence of an additional X-ray field, would be 
useful for a more quantitative evaluation of this scenario. 



6.3. Shock Heating 

The HE CO component is fit by a C-shock model 
with preshock density — 10^ cm~'^, veloc- 

ity Vs = 30 kms""'^, and cross-section A ~ (21 
pc)2 (|Kaufman fc Neufeldl ir996) (Figure (H}:). As dis- 
cussed in sections 14.21 and 15.31 jet-ISM interactions in 
the H2 ring are unlikely to produce the FIR CO emis- 
sion, but a jet-induced shock in the northern streamer is 
a more plausible source of the HE CO lines. The radio jet 
changes direction in the vicinity of the brightest clump in 
the northern streamer, evidence that the fast moving, low 
density material in the jet is collidi ng with and being di- 
verte d by the dense molecular gas (jMuller Sanchez et al.l 
|2009( ). This interaction should produce shocks in the 
molecular material, and the pres ence of H2O masers i n 
this cloud supports this picture (jGallimore et al.ll2004D . 
The modeled A ^ {21 pc)^ cross-section matches the ~14 
pc size of this clump, consistent with this scenario. The 
je t mechanical power estima ted in the bow shock model 
of IWilson fc UlvestadI ([l987h is ~2 x 10^ Lq, a factor of 
^3 larger than the mechanical luminosity L = l/2pv^A 
of the shock model discussed here. Jet-induced shocks 
are therefore energetically feasible, although would re- 
quire an efficient mechanism for converting kinetic en- 
ergy at high velocities into slow molecular shocks. Al- 
ternatively, shocks in the bluest regions of the H2 ring, 
possibly associated with the ring expansion, could also 
generate the HE CO. 

An important constraint on this shock modeling is 
the need to match the CO brightness while not over- 
producing the ro-vibrational H2 emission. For the 
C-shock model discussed above the predicted H2 1- 
S(l) flux is a factor of ~70 larger than ob served 
in the northern streamer (jMuller Sanchez et al.l [2009) , 
and ~7 larger than in the bright eastern clump in the 
H2 ring (TG alliano & AUoin 2002). In this region of 
the iKaufman fc Neufeldl (|1996f ) model parameter space 
the H2/CO intensity ratio is a steep function of shock 
velocity, while the CO emission is sensitive to both the 
velocity and preshock density. Finely-tuning Vg and uq 
(going to lower Vs and higher no) may therefore yield 
a solution that reduces the H2 emission while conserv- 
ing the CO line SED. Extinction of the 2 /im H2 emis- 
sion may a lso be important, part i cularl y in the northern 
streamer. iMiiller sknchez et al.l ()2009[ ) estimate a col- 
umn of A^H ~ 8 X 10^'* cm~^ in the southern streamer. 
Assuming a similar value for the northern streamer, and 
for Ak/Nh = (2 X 10^2 cm-2)-\ the resulting Ak ~ 400 
is more than sufficient to provide the necessary atten- 
uation in a mixed dust model. A J-shock model for 
the HE CO would also generate much weaker H2 emis- 
sion ( Hollenbach fc McKeel ll989'). although the required 
cross-section would then be a factor of ~5 higher than for 
a C-shock, and no longer match the northern streamer 
size. In sum, a number of scenarios may be envisioned to 
attribute the HE CO lines to a shock in the H2 ring or 
the northern streamer, while also satisfying the H2 1-0 
S(l) brightness. 

6.4. Summary and Discussion 

We conclude that the HE component most plausibly 
arises from X-ray or shock heated gas in the northern 
streamer, or shock heated gas in the H2 ring. For the 
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ME component, we argued in section [53] that energetic 
considerations and the ancihary evidence for an XDR 
chemistry favored X-ray over shock heating, but the sit- 
uation for the HE component is less clear. The amount 
of mechanical power available in the jet or other sources 
is indeed uncertain, but the strong evidence for a jet 
interaction with the northern streamer suggests that jet- 
induced shocks in this cloud be given full consideration. 
Additionally, the 0H+ and H2O+ emission lines that 
pinpoint an XDR chemistry are emitted close to the 
galaxy systemic velocity, with no detected emission at 
the blueshifted velocity of the HE CO lines (Figure |6]). 
Thus while nuclear X-rays may dominate the energet- 
ics and chemistry in the CND, an origin of the HE CO 
lines in a small amount of shock-heated gas must also be 
considered. 



7. NON-DETECTIONS, AND CONSTRAINTS ON THE 
NUCLEAR OBSCURATION 

In addition to the detected Jupper = 14—24 transitions, 
our upper limits to the Jupper < 50 lines provide a use- 
ful constraint on the nuclear molecular component. One 
of the earliest discussions o f high-J CO in extragalac- 
tic sources was presented in lKrolik fc LeppI ([1989), who 
modeled the emission expected from the Compton-thick 
material obscuring type 2 AGN. These authors calcu- 
lated the CO emission from an A^h = 10^"* cm~^ cloud 
located ~1 pc from a luminous {Lx ~ 10'*'' ergs"') 
hard X-ray source, in which the heating is dominated 
by Compton scattering of 10 — 100 keV photons. In this 



model the CO SED scales as Jupper up to Jupper ~ 58, 
and the absolute line luminosities are proportional to 
the total absorbed 10 — 100 keV luminosity {fahsLx44', 
in units of lO''^ ergs"'). Several of our upper limits 
to CO transitions with Ju 



pper 



30 — 50 independently 
place an upper limit to this component corresponding to 
fahsLxAA ^ 0.1, while the non-detection of the CO(44-43) 
line achieves the most stringent limit of fahsLxM ^ 0.09 
(Figure [8^). This limit may be reduced by stacking 
the individual non-detections. We have obtained such 
a stack by first scaling the spectrum of each undetected 
line by (44/ Jupper)^, which references each spectrum to 
that of 00(44-43) under the assumption that the CO 
fluxes scale as Jupper- We then calculated a weighted 
average of 8 lines with low noise that fall in clean spec- 
tral regions, and binned to 600 kms~'. The result is 
shown in Figure [TOl Th e upper limit to this s tacked line 
pushes the limit of the iKrolik fc LeppI ()1989D model to 
/absi,44 < 0.038 (Figure |3i). 

Hard X-ray surveys co n ducted with INTE- 
GRAL /IBIS (lMa,lizia et al.l [20091 ) and Swift- 
BAT (jBurlon et al.l 120 lit) indicate -20 - 25% of 
AGN are Compton thick, and in the spirit of uni- 
fication we adopt this as the covering factor of the 
Comp ton-th ick material in NGC 1068. Ilwasawa et al.l 
(|1997| ) and iColbert et cH (|2002[ ) model the reflected 
X-ray spectrum of NGC 1068, and estimate intrinsic 
luminosities of L2-iokeV — lO^'^"^'^-^ ergs"' (corrected 
to the d = 14.4 Mpc used here). A comparison of 
the measured [OIII] A5007 line luminosity with a 
log([OIII]/Lx) = —1.89 ratio established for type 1 
Seyferts yields i2-iokcV = lO'''^'^ ergs"', while a com- 
parison of the estimated Lboi with L2-iokev/-^boi ~ 0.1 
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Fig. 10. — Baseline-subtracted stack of 8 transition s with 
Jupper = 34 — 47, for comparison with the IKrolik fc Lepd 1)19891 ) 
models. Dashed lines mark the ±lcr uncertainty. 



yields i2-iokeV ^ 

204 3.6 ergs-' (IColbert et al.l [2002L 
and references therin). iMelendez et al.l ( 20081) derive a 
relation between i[oiv]26Atm and L2-iokeV for a sample 
of hard (14-195 keV) X-ray selected type 1 Seyferts. 
Using this relation, along with i[oiv]26/im = 4.7 x 10"" 
ergs"' ([Sturm et al.l I2002f) . yields L2-iokeV ~ 10''^ 
ergs"'. Given this range of estimates we adopt 
i2-iokeV = lO''^'^ ergs"', and increase this by 1.43 
to correct to the 10 — 100 keV range (for cx i^~'). 
Assuming 25 % of this powe r is ab sorbed in the phase 
modeled by IKrohk fc Lerol ([1989D gives an expected 
fiihsLx44 = 0.11, a modest factor of ^-^2.9 higher than our 
measured upper limit. Given the uncertainties in the 
AGN luminosity and covering factor of the Comption- 
thick material in NGC 1068, as well as uncertainties 
in the model itself, this factor of ~2.9 overprediction 
is too small to rej ect the model of a parsec-scale torus 
envisioned bv iKrolik fc Lepd (|1989D . This non-detection 
may, however, provide a useful constraint on the more 
recent set of detailed tor us models (e.g., iHonig et al.l 
l200l[Nenkova et al.ll2008l) . 

8. SUMMARY AND FUTURE WORK 

We have detected 10 CO transitions in the Jupper = 
14 - 24 range from the central 10" (700 pc) of NGC 
1068, and obtained sensitive upper limits to most other 
transitions up to Jupper < 50. These are the first ex- 
tragalactic detections of FIR CO, which represent a new 
probe of excited molecular gas in Seyfert nuclei. The 
detected transitions are modeled as arising from 2 differ- 
ent components: a moderate excitation (ME) component 
close to the galaxy systemic velocity, and a high excita- 
tion (HE) component that is blueshifted by —70 kms~'. 
Our main results are as follows: 

1) Using a two component LVG model we derive nH2 ^ 
10^-^ cm-^, Tkin 146 K, and Mh2 lO^'^ Mq for the 
ME component, and 71^2 ~ lO'^'^ cm"^, ~ 435 K, 
and Mh2 ^ lO^'* Mq for the HE component. The la un- 
certainties on the derived temperatures are '-^±0.25 dex, 
while for density and mass this is — ±1 dex. Extend- 
ing the measured CO line SED to lower- and higher- J 
lines would likely reduce these uncertainties, as would a 
joint analysis with H2, OH, H2O, and other complemen- 
tary molecular tracers. We will follow both approaches 
in forthcoming papers. 
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2) These two components are denser than the gas 
traced with milhmeter CO observations, and the HE (and 
possibly the ME) component is also warmer. The ME 
component makes a non-negligible contribution to the 
nuclear mass budget, although large uncertainties in the 
masses estimated from both the FIR CO lines and from 
CO(l-O) prevent a more quantitative statement. 

3) Comparing the FIR CO line profiles with those 
of high spatial and spectral resolution observations of 
C0(2-l) and H2 1-0 S(l) allows a first estimate of the 
origins of the ME and HE components within the central 
10". Good matches are found with H2 1-0 S(l), which 
for the ME component suggests an origin in the ~200 
pc diameter H2 ring. The blueshifted HE lines may also 
be consistent with an origin in the bluest regions of the 
H2 ring, but are better matched to the clump of infalling 
molecular gas ^40 pc north of the AGN. 

4) The ME component is nicely consistent with models 
of X-ray heating of gas in the CND. A shock model is 
also possible, although due to the uncertainties in the 
amount of mechanical power available for dissipation in 
slow shocks, and the evidence for X-ray driven chemistry 
in the CND, we favor the X-ray heating scenario. Far-UV 
heating is unlikely. 

5) The HE component is also consistent with either 
X-ray or shock heating. X-ray heating would best fit 
with an origin in the cloud ~40 pc north of the AGN, 
supporting the results of the line profile matching (point 
3). Shocks triggered by the interaction of the radio jet 



with this clump, or arising from the H2 ring, are also 
plausible. Far-UV heating is unlikely. 

6) Our non-detections of Juppor = 30 — 50 lines place a 
constraint on any gas associated with a putative parsec - 
scale torus. Our upper limit to the lKrolik fc Leppll) 19891 ) 
prediction is a factor of ^-^2.9 lower than the expected sig- 
nal, although due to the uncertainties involved in apply- 
ing this model, this non-detection is insufficient to rule 
out the torus paradigm. The inclusion of a gas phase in 
the current set of clumpy torus models, and a compar- 
ison of the predicted CO line SED with our detections 
and upper limits, would be a useful next step. 
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